| Name | LAI (-) | Canopy Coverage (-) |
|---|---|---|
| Sparse | 1.59 | 0.73 |
| Mixed | 1.86 | 0.78 |
| Closed | 2.11 | 0.82 |
Snow Interception Relationships with Meteorology and Canopy Structure in a Subalpine Forest
Authors:
A. Cebulski1 (ORCID ID - 0000-0001-7910-5056)
J.W. Pomeroy1 (ORCID ID - 0000-0002-4782-7457)
1Centre for Hydrology, University of Saskatchewan, Canmore, Canada
Corresponding Author: A. Cebulski, alexcebulski@gmail.com
Abstract: Subcanopy snow accumulation models differ in how snow interception and ablation processes are represented and thus their application to diverse climates and forest types is uncertain. Existing parameterizations of initial snow interception before unloading include inherently coupled canopy snow accumulation and ablation processes. This leads to difficulty in diagnosing processes and adding possible errors to simulations when incorporated as canopy interception routines in models that already account for canopy snow ablation. This study evaluates the theory underpinning these parameterizations using in-situ meteorological data, high-temporal resolution point-scale throughfall measurements, and fine-scale aerial lidar measurements of throughfall and canopy metrics collected from two subalpine forest plots in the Canadian Rockies. Contrary to existing theories, no association of canopy snow load or air temperature with interception efficiency was observed. Instead, forest structure emerged as the primary factor governing snow accumulation. A wind-driven snowfall event demonstrated that non-vertical hydrometeor trajectories can significantly increase interception. Prediction of interception efficiency for this event improved dramatically when adjusted for hydrometeor trajectory angle based on a wind speed at one-third of the canopy height. Snow-leaf contact area showed a high sensitivity to wind speed, increasing by up to 95% with a 1 m s-1 wind speed. The study proposes a new parameterization which calculates snow interception efficiency as a function of snow-leaf contact area adjusted for hydrometeor trajectory angle. This new parameterization successfully estimated subcanopy snow accumulation for a snowfall event at two forest plots measured using lidar and snow surveys. By separating canopy snow ablation from snow interception processes, this new model offers potentially improved prediction of subcanopy snow accumulation when combined with canopy snow ablation parameterizations.
Keywords: snow interception, throughfall, ablation, forest, snowpack, lidar, process-based modelling
Introduction
Over half of North America’s snow-covered zone is covered by forests (Kim et al., 2017), significantly impacting the accumulation and redistribution of snowpacks and subsequent snowmelt runoff. Essery et al. (2003) estimated that 25–45% of annual snowfall may be lost to the atmosphere due to sublimation of snow intercepted in forest canopies globally. Snow intercepted in the canopy can sublimate and melt at much higher rates than the subcanopy snowpack (Floyd, 2012; Lundberg & Hallidin, 1994; Pomeroy et al., 1998), reducing the amount of snow available for runoff. Forest thinning efforts aimed at limiting sublimation losses to increase snowmelt runoff do not always lead to a corresponding increase in spring streamflow (Golding & Swanson, 1978; Harpold et al., 2020; Pomeroy et al., 2012; Troendle, 1983). This may be due to increased ablation rates when forest cover is reduced, desynchronization of snowmelt, and sub-surface hydrology interactions (Ellis et al., 2013; Musselman et al., 2015; Pomeroy et al., 1997; Safa et al., 2021; Varhola et al., 2010). Vegetation structure controls the partitioning of snowfall into throughfall and interception, and thus governs the quantity of snow subject to sublimation from the canopy (Hedstrom & Pomeroy, 1998; Storck et al., 2002). Due to the significant impact of forest cover on snow accumulation and ablation, and sparse or absent monitoring networks for subcanopy snow accumulation (Rittger et al., 2020; Vionnet et al., 2021), land management, ecological conservation and water resource decisions rely on robust models of snow redistribution to estimate past, current and future subcanopy snowpacks.
Hedstrom & Pomeroy (1998), working in the cold continental boreal forest, proposed that initial snow interception efficiency was controlled by the maximum canopy load which itself was a function of leaf area index and air temperature. Unloading was found to be an exponential function of time and observed only days or weeks after the interception event. Storck et al. (2002), working in temperate coastal forests, emphasized the role of air temperature in controlling the maximum canopy snow load. Gelfan et al. (2004) demonstrated accurate subcanopy snowpack simulations at study sites in Russia by treating the Hedstrom & Pomeroy (1998) and Storck et al. (2002) parameterizations separately while using a step-based function to choose either parameterization based on temperature. A similar parameterization in the Cold Regions Hydrological Model (Pomeroy et al., 2022) has shown strong performance at sites across Canada, northern United States, Switzerland, and Spain. However, overestimation of subcanopy snow accumulation was reported by Lundquist et al. (2021) and Lumbrazo et al. (2022) when combining the Hedstrom & Pomeroy (1998) routine with ablation parameterizations from different studies (e.g., Roesch et al., 2001). The coupling of ablation processes within existing snow interception models (Hedstrom & Pomeroy, 1998; Storck et al., 2002) may contribute to the over representation of ablation processes when combined with other canopy snow ablation parameterizations. Additional observations of snow interception that exclude ablation processes could help determine the applicability of the interception theories proposed by Hedstrom & Pomeroy (1998) and Storck et al. (2002). Hedstrom & Pomeroy’s (1998) theory also suggests that moderate wind speeds, which can result in more horizontal hydrometeor trajectories and increase the snow-leaf contact area and interception efficiency at the plot scale. This association has also been shown in rainfall interception studies (i.e., Herwitz & Slye, 1995; Van Stan et al., 2011) to decrease throughfall of rain. Despite this importance for rainfall, the relationship proposed by Hedstrom & Pomeroy (1998), has typically not been included in snow accumulation models (Clark et al., 2020; Mahat & Tarboton, 2014) and empirical testing of this relationship is lacking.
The objective of this paper is to evaluate the theories underlying existing snow interception models using high spatial and temporal resolution measurements of subcanopy snow accumulation for events with minimal canopy snow ablation. These new observations are investigated to address the following research questions:
Are the existing theories regarding the relationships between meteorology and forest structure and snow interception supported by in-situ observations?
Is snow interception influenced by non-vertical hydrometeor trajectory angles over a wind-driven snowfall event?
To what extent can these findings inform the development of a new parameterization for snow interception?
Theory
Snow Interception
During calm snowfall periods, where canopy snow wind redistribution processes and ablation processes such as sublimation and melt/drip can be assumed negligible, the canopy snow load, \(L\) (kg m-2) can be estimated from the mass balance:
\[ \frac{dL}{dt} = q_{sf} - q_{tf} - q_{unld}(L) - q_{drip}(L) - q_{wind}^{veg}(L) - q_{sub}^{veg}(L) \tag{1}\]
where \(q_{sf}\) is the snowfall rate (kg m-2 s-1), \(q_{tf}\) is the throughfall rate (kg m-2 s-1), \(q_{unld}\) is the canopy snow unloading rate (kg m-2 s-1), \(q_{drip}\) is the canopy snow drip rate due to canopy snowmelt (kg m-2 s-1), \(q_{wind}^{veg}\) is the wind transport rate in our out of the control volume (kg m-2 s-1), and \(q_{sub}^{veg}\) is the intercepted snow sublimation rate (kg m-2 s-1).
During periods with low air temperatures and low wind speeds where \(q_{unld}\), \(q_{drip}\), \(q_{wind}^{veg}\), and \(q_{sub}^{veg}\) can be assumed negligible.
Interception efficiency, \(\frac{I}{P}\) (-), which is the fraction of snowfall intercepted over \(\Delta t\) can then be calculated as:
\[ \frac{I}{P} = \frac{\frac{dL}{dt}}{q_{sf}} \tag{2}\]
and throughfall, \(q_{tf}\) can be calculated as:
\[ q_{tf} = \left(1 - \frac{I}{P} \right) \cdot q_{sf} \tag{3}\]
Hydrometeor Trajectory Angle
The trajectory angle, \(\theta_h\) of a hydrometeor as the departure in degrees (°) from a vertical plane (i.e., 0° for vertical snowfall), is shown in Herwitz & Slye (1995) to be calculated as:
\[ \theta_h = \arctan \left(\frac{x_h(u_z)}{v_h(D_h)}\right)*\frac{180}{\pi} \tag{4}\]
where \(v_h(D_h)\) is the terminal fall velocity of the hydrometeor (m s-1), which is a function of the hydrometeor diameter, \(D_h\) and \(x_h(u_z)\) is the horizontal velocity of the hydrometeor (m s-1) which is a function of the within canopy wind speed, \(u_z\) at height above ground, \(z\), if the hydrometeors are assumed to be following fluid points in the atmosphere.
Data and Methods
Study Site
This study was conducted at Fortress Mountain Research Basin (FMRB), Alberta, Canada, -115° W, 51° N, a continental headwater basin in the Canadian Rockies (Figure 1). Data from this study was collected between October 2021 and July 2023 within and surrounding two forest plots adjacent to the FMRB Powerline Station (PWL) and Forest Tower Station (FT) at ~2100 m above sea level as shown in Figure 1. The average annual precipitation at PWL Station from 2013 to 2023 was 1045 mm, with the peak annual snow water equivalent (SWE) reaching 465 kg m-2, typically occurring in late April. The PWL and FT forest plots include discontinuous stands of 70% subalpine fir (Abies lasiocarpa) and 30% Engelmann spruce (Picea engelmannii) (Langs et al., 2020). The PWL plot is located 120 m to the northwest of FT station and contains a forest clearing with a diameter of ~12 m, surrounded by a closed canopy. The canopy coverages of the two forest plots are 0.51 and 0.29 and the winter leaf area indices are 2.07 and 1.66 for PWL and FT respectively. The average height of the canopy surrounding the plot to the east of the PWL station is 10.5 m and surrounding the forest plot around the FT Station is 7.1 m. The forest of the FT plot has a discontinuous canopy without artificial clearings. In August of 1936, the majority of vegetation in FMRB burned during a large forest fire that affected most of the Kananaskis Valley (Fryer et al., 1988). Following the fire, the forest within the PWL and FT forest plots has naturally regenerated, though some trees have been removed for road clearing and creation of a snow study plot.
Meteorological Measurements
Measurements of air temperature and relative humidity (Vaisala model HMP155A), wind speed and direction (RM Young model 86000 2-D ultrasonic anemometer) were made 4.3 m above the ground at FT station (Figure 1). Wind speed measurements from a 3-cup anemometer (Met One model 014A), installed adjacent to the 2-D ultrasonic anemometer at 4.3 m, were used for gap filling wind speed. Additional wind speed measurements were collected by two 3-D sonic anemometers (Campbell Scientific CSAT3) installed at 2 m (raised to 3 m February 2022) and 13.5 m above the ground at FT station. Average wind speeds at these four heights, for periods where the instruments were known to be clean of snow, were found to follow a logarithmic relationship and a wind profile was fitted using the Prandtl-von Kármán log-linear relationship:
\[ \overline{u} = \frac{u_*}{k} ln(\frac{z - d_0}{z_0}) \tag{5}\]
where \(\overline{u}\) is average wind speed (m s-1) at height, \(z\) (m) above the ground, \(u_*\) is the friction velocity (m s-1), \(d_0\) is the displacement height (m), \(z_0\) is the roughness length of momentum (m), and \(k\) is the dimensionless von Kármán Constant (0.4).
To determine the displacement height and roughness length parameters the function “optim” from the stats R package (R Core Team, 2024) was used. The parameters for the wind speed profile include a roughness length of 0.50 m, displacement height of 0.58 m. At PWL station, the snowfall rate was measured by an Alter-shielded OTT Pluvio weighing precipitation gauge 2.6 m above ground, corrected for undercatch following phase correction by Harder & Pomeroy (2013) and catch efficiency by Smith (2007). Wind speed for undercatch correction was measured by a 3-cup anemometer (Met One model 014A) at a height of 2.6 m at PWL station. An optical disdrometer (OTT Parsivel2) provided measurements of hydrometeor particle size and vertical velocity. All measurements were recorded at 15-min intervals using Campbell Scientific dataloggers, except the Parsivel2 which was recorded at 1-minute intervals by an onsite computer.
Lysimeter Measurements
Three subcanopy lysimeters (SCLs) were installed surrounding the FT Station (Figure 1) to provide 15-minute interval measurements of sub-canopy snowfall and canopy snow load as in MacDonald (2010). Figure 2 shows the three SCLs which consisted of a plastic horse-watering trough with an opening of 0.9 m2 and depth of 20 cm suspended from a load cell (Intertechnology 9363-D3-75-20T1) attached to an aluminum pipe connected between two trees. For 26 distinct snowfall events, where canopy snow ablation rates were deemed negligible, the throughfall rate, \(q_{tf}\), was calculated by dividing the load cell weight by the cross-sectional area of the SCL opening and calculating the rate of change at 15-minute intervals. Canopy snow load and interception efficiency was quantified at the same 15-minute intervals during these events using ?@eq-dwdt-ode and Equation 2, incorporating measurements of \(q_{tf}\) from the SCLs and \(q_{sf}\) from the PWL snowfall gauge. Timelapse imagery, mass change on a weighed tree lysimeter “hanging tree” (Pomeroy & Schmidt, 1993) and in-situ observations were used to ensure the ablation of snow intercepted in the canopy or snow on the ground was minimal over each of the selected snowfall events. Additionally, the \(q_{tf}\) measurements were filtered to include observations with a snowfall rate > 0 kg m-2 hr-1, throughfall rate > 0.05 kg m-2 hr-1 and a snowfall rate greater than the subcanopy lysimeter throughfall rate to minimize observations with unloading. The weighed tree lysimeter, a live subalpine fir (Abies lasiocarpa) tree suspended from a load cell (Artech S-Type 20210-100) measured the weight of canopy snow load, \(L_{wt}\) (kg). The weight of snow intercepted on the weighed tree was scaled to an areal estimate of canopy snow load (\(L\), kg m-2) using measurements of areal throughfall (kg m-2) from manual snow surveys and snowfall from the PWL Station snowfall gauge (see description of method in Pomeroy & Schmidt, 1993). The canopy structure surrounding three SCLs is shown in Figure 2 and was measured using hemispherical photography (Nikon Coolpix 4500 and EC-F8 hemispherical lens) and the hemispheR R package Chianucci & Macek (2023). The leaf area index and canopy coverage from hemispherical photo analysis is shown in Table 1.
UAV-Lidar Data Collection Processing
Two uncrewed aerial vehicle (UAV) lidar surveys were conducted before and after a 24-hour snowfall event that occurred between March 13th and March 14th, 2023 to facilitate the measurement of snow accumulation and canopy structure metrics. The UAV (FreeFly Alta X) payload included a REIGL miniVUX-2 airborne laser scanner, an Applanix APX-20 inertial measurement unit (IMU) and global navigation satellite system (GNSS). The UAV was flown 90 m above the ground at a speed of 3 m s-1 following the path shown in Figure 1. A detailed description of the UAV, payload, and flight settings is provided in the supporting information. The methods outlined by Harder et al. (2020) and Staines & Pomeroy (2023) were incorporated to reconcile survey lidar, IMU and GNSS data. A vertical offset of up to 6 cm between UAV-lidar flight lines was observed in the resulting point clouds on March 13th and 14th, 2024 and was attributed to IMU position drift. This offset between flight lines was corrected using the BayesStripAlign software v2.24 (BayesMap Solutions, 2024). After strip alignment, the mean elevation bias (lidar minus GCP) was 0.000 m and the RMS error declined from 0.055 m to 0.038 m March 13th and changed from 0.033 m to 0.029 m on March 14th. The point cloud density ranged from ~1200 returns m2 in open clearings to ~2200 m2 in sparse forest for both the March 13th and 14th surveys after all flight paths were combined. Quality control, ground classification and calculation of the change in between two UAV-lidar point clouds was conducted using the LAStools software package (LAStools, 2024). More details on the UAV-lidar processing workflow are provided in the supporting information.
Snow Surveys
In-situ Snow Depth and Density
In-situ fresh snow surveys provided measurements of subcanopy throughfall depth and density following the transects shown in Figure 1. Twelve fresh snow surveys (six pre- and post-snowfall event pairs) which had minimal ablation and redistribution between pre and post surveys at 30 locations were selected to upscale the weighed tree snow load. When conditions allowed for a UAV-lidar flight, the in-situ snow surveys were conducted following the UAV-lidar flight to assess the accuracy of the throughfall measurements and provide a fresh snow density for the calculation of SWE (kg m-2). A 1000 cm3 Perla snow density wedge sampler (RIP Cutter, https://snowmetrics.com/shop/rip-1-cutter-1000-cc/) was used to measure the density of the fresh snow layer, \(\overline{\rho_{tf}}\) (kg m-3) from snow pits. Throughfall depth measurements, \(\Delta HS\) were converted to SWE using the following equation:
\[ \Delta SWE_{tf} = \Delta HS \cdot \overline{\rho_{tf}} \tag{6}\]
Differential GNSS rover coordinates, with ± 2.5 cm 3D uncertainty, were taken at each snow sampling location so the locations could be queried later from the UAV-lidar rasters. If a pre-event crust layer was present, the depth of post event fresh snow accumulation above the crust layer was interpreted as throughfall over the event. In the absence of a defined crust layer, the difference in pre- and post-event snow depth to ground was interpreted as event throughfall.
UAV-Lidar Snow Depth
Two UAV-lidar surveys, one prior to a snowfall event on March 13, 2023 at 10:00 CST and another following snowfall on March 14, 2023 at 11:00 CST enabled fine-scale analysis of snow accumulation and canopy structure within the FT and PWL forest plots. This period was selected based on two criteria: 1) it provided sufficient cumulative snowfall to result in a low relative error in UAV-LiDAR measured throughfall, and (2) minimal redistribution and ablation was observed, as confirmed by the SCLs, weighed tree, and time-lapse imagery. The change in elevation between the two UAV-lidar surveys was interpreted as the increase in snow accumulation, \(\Delta HS\) over the snowfall event. Further details on the generation of 25 cm horizontal resolution \(\Delta HS\) rasters from the UAV-lidar point clouds are provided in the supporting information.
UAV-Lidar Canopy Metrics
To characterize the canopy structure, the voxel ray sampling (VoxRS) methodology for lidar data analysis was employed, as developed by Staines & Pomeroy (2023), for both UAV-lidar surveys. This method was chosen for its ability to provide canopy metrics that are less sensitive to the inherent non-uniform nature of lidar sampling data, which often results from beam occlusion in vegetation and leads to reduced points near the ground. The VoxRS algorithm is publicly available at https://github.com/jstaines/VoxRS. The canopy products produced from VoxRS here include canopy contact number, the mean theoretical number of canopy contacts for a given ray, and radiation transmittance (\(\tau\)), all dimensionless. See supporting information in Staines & Pomeroy (2023) for details on how these metrics are computed. The fraction of snow-leaf contact area per unit area of ground proposed by Hedstrom & Pomeroy (1998), and hereafter called leaf contact area (\(C_p\)), was calculated as:
\[ C_p(C_c, \theta_h, L) = 1-\tau \tag{7}\]
\[ C_p(C_c, \theta_h, L) = \begin{cases} 1-\tau,& \text{if } \theta_h> 0°\\ 1-\tau \approx C_c , & \theta_h= 0° \end{cases} \tag{8}\]
where \(C_p\) is a function of the canopy coverage \(C_c\), \(\theta_h\) and \(L\). \(C_p\) is approximately equal to canopy coverage (\(C_c\)) for vertical snowfall trajectories. However, for non-vertical snowfall \(1 > C_p > C_c\).
Statistics and Regression Models
To determine how forest structure was associated with interception efficiency at different azimuth and zenith angles over the March 13-14 snowfall event, each portion of the hemisphere at each grid location was considered. The relationship between interception efficiency and canopy contact number was found to be linear and thus the Pearson Correlation Coefficient, \(\rho_p\) was calculated using the ‘stats’ package in R (R Core Team, 2024) to quantify the association between a single raster of interception efficiency and the 32,760 rasters containing the canopy contact number hemisphere for each portion of the hemisphere (azimuth [0°, 1°, …, 359°], zenith angle [0°, 1°, …, 90°]) for each of the 25 cm grid cells across the FT and PWL forest plots.
Linear and non-linear regression models were developed to assess relationships in the observed data. Linear models were fitted using ordinary least squares regression via the ‘lm’ function from the R ‘stats’ package (R Core Team, 2024) to analyze two relationships: (1) between interception efficiency and meteorological variables and (2) between interception efficiency and leaf contact area. The latter was forced through the origin based on the theoretical justification that the dependent variable should be zero when the independent variable is zero. Kozak & Kozak (1995) noted, the default R2 value provided for least squares models forced through the origin by many statistical packages can be misleading. Therefore, these R2 values were adjusted using Equation 10 in Kozak & Kozak (1995). Non-linear models were fitted using non-linear least squares (nls) regression via the ‘nls’ function in ‘stats’ package in R.
Results
The influence of meteorology on snow interception
Canopy snow load was estimated for 26 snowfall events, and increased linearly with cumulative event snowfall without evidence of reaching a maximum (Figure 3). The air temperature ranged from -24.5°C to 1°C (Table 2), wind speeds at 4.3 m height ranged from calm to 4.6 m s-1, and wind direction was predominately from the southwest during snowfall (Figure 4). Missing canopy snow load measurements in Figure 3 for certain troughs during specific events was caused by damage to the subcanopy lysimeter wiring due to animals and heavy snow loads.
The mean event air temperature was weakly negatively associated (p < 0.05) with mean event interception efficiency at the mixed canopy, but not associated at the closed and sparse canopies (Figure 5). Cumulative event snowfall was not associated with mean event interception efficiency (p > 0.05). Event mean wind speed was positively associated with interception efficiency for the sparse (R2 = 0.1, p > 0.05) and closed (R2 = 0.2, p < 0.05) canopies, both with limited canopy openings (Figure 2, A & C) towards the prevailing wind direction (Figure 4). The mixed canopy was not associated with wind speed (p > 0.05).
| Dependent Variable | SCL Name | \(C_c\) | Adjusted \(R^2\) | \(p\)-value | \(n\) |
|---|---|---|---|---|---|
| Cumulative Snowfall (kg m⁻²) | Mixed | 0.78 | 0.03 | 0.20 | 26 |
| Air Temperature (°C) | Mixed | 0.78 | 0.14 | 0.03 | 26 |
| Wind Speed (m/s) | Mixed | 0.78 | 0.01 | 0.27 | 26 |
| Cumulative Snowfall (kg m⁻²) | Closed | 0.82 | -0.05 | 0.73 | 20 |
| Air Temperature (°C) | Closed | 0.82 | 0.01 | 0.30 | 20 |
| Wind Speed (m/s) | Closed | 0.82 | 0.19 | 0.03 | 20 |
| Cumulative Snowfall (kg m⁻²) | Sparse | 0.73 | -0.04 | 0.57 | 19 |
| Air Temperature (°C) | Sparse | 0.73 | -0.03 | 0.52 | 19 |
| Wind Speed (m/s) | Sparse | 0.73 | 0.11 | 0.09 | 19 |
Fifteen minute interval measurements of mean interception efficiency and air temperature were not associated, despite significant relationships (R2 < 0.03, p < 0.05) for the sparse and mixed canopies, due to low predictive power (Figure 6, A). The average interception efficiency across differing bins of air temperature also do not show any systematic trend (Figure 6, A). However, a significantly greater median interception efficiency (p < 0.05) was found for binned measurements with air temperatures below -6 °C compared to those with warmer air temperatures using non-parametric Wilcoxon signed rank test.
No association (R2 < 0.09, p < 0.05 for all three canopies) was found between wind speed measured at the FT Station and interception efficiency for the 15-minute measurements (Figure 6, B). The binned data show an increasing trend in interception efficiency with increasing wind speed for the sparse and closed canopies (Figure 6, B). A comparison of interception efficiencies binned for low (< 1 m s-1) and high (> 1 m s-1) wind speeds by the Wilcoxon signed rank test, showed that high wind speeds had significantly higher (p < 0.05) median interception efficiencies compared to the low wind speed bins for the closed and sparse canopy. Conversely, the Wilcoxon test showed the mixed canopy, which had an opening in the canopy towards the prevailing wind direction (Figure 2, B), had significantly higher (p < 0.05) median interception efficiencies for the low wind speed bins.
Interception efficiency measured at 15-minute intervals showed no association with the initial canopy load (Figure 6, C). The binned data show a small increase in interception efficiency with initial snow load for snow loads less than 7 kg m-2 for all canopies, but interception efficiency declined for initial snow loads greater than 7 kg m-2 at all canopies. Though this was inconsistent for the mixed canopy. A comparison of snow interception for low (< 10 kg m-2) and high (> 10 kg m-2) initial canopy snow load bins using the Wilcoxon rank-test showed the low canopy snow loads had significantly greater (p < 0.05) median interception efficiencies than those with high initial canopy snow loads.
| Dependent Variable | SCL Name | \(C_c\) | Adjusted \(R^2\) | \(p\)-value | \(n\) |
|---|---|---|---|---|---|
| Air Temperature (°C) | Mixed | 0.78 | 0.03 | 0.00 | 985 |
| Air Temperature (°C) | Closed | 0.82 | 0.00 | 0.07 | 618 |
| Air Temperature (°C) | Sparse | 0.73 | 0.01 | 0.02 | 603 |
| Wind Speed (m/s) | Mixed | 0.78 | 0.02 | 0.00 | 985 |
| Wind Speed (m/s) | Closed | 0.82 | 0.04 | 0.00 | 618 |
| Wind Speed (m/s) | Sparse | 0.73 | 0.09 | 0.00 | 603 |
| Initial Canopy Snow Load (kg m⁻²) | Mixed | 0.78 | 0.00 | 0.45 | 972 |
| Initial Canopy Snow Load (kg m⁻²) | Closed | 0.82 | 0.05 | 0.00 | 607 |
| Initial Canopy Snow Load (kg m⁻²) | Sparse | 0.73 | 0.03 | 0.00 | 592 |
The influence of forest structure on snow interception
UAV-lidar measurements of throughfall and canopy structure metrics provide insights on how the forest canopy influenced subcanopy snow accumulation during a wind-driven snowfall event between March 13th and 14th 2023. This event totaled 28.7 kg m-2 of snowfall at PWL station and was characterized by a transition from low rates of snowfall and air temperatures near 0°C to higher rates of snowfall by late afternoon on March 13th coinciding with air temperatures around -2.5 °C. An average wind speed of 1.3 m s-1 and direction of 188° was observed 4.3 m above the ground at FT Station. The log-linear wind speed profile shown in Figure 7 provided a good fit to observed wind speeds at 2, 3, 4.3 and 13.5 m above the ground and shows Cionco’s (1965) exponential function was not appropriate for this sparse canopy. The heavy snowfall over this event covered the two eddy covariance systems at FT station with snow, and thus this wind profile was only assessed using the wind speed measured at 4.3 m for this event (Figure 7). Figure 7 shows predicted hydrometeor trajectory angles at varying heights, calculated using Equation 4 and the mean observed hydrometeor terminal velocity observed over the event of 0.9 m s-1. An average wind speed of 1.6 m s-1 and direction of 188° was calculated by integrating the wind speed from the surface to the mean canopy height of FT plot. The corresponding trajectory angle calculated using Equation 4 from this integrated wind speed was 61.5°.
Throughfall depth measurements by UAV-lidar aligned well with 28 in-situ manual measurements with a mean bias of -0.001 m and RMSE of 0.024 m. More details on the accuracy of UAV-lidar snowdepth measurements are provided in the supporting information section. Figure 8 shows the spatial distribution of throughfall and interception efficiency at the PWL and FT forest plots. Reduced throughfall and greater interception efficiency was observed on the north (lee) side of individual trees, which may be due to non-vertical hydrometeor trajectories caused by the steady southerly winds observed over this event. Visual observations on March 13th and 14th confirmed non-vertical hydrometeor trajectories and increased canopy snow loads were observed on the windward side of individual trees. This effect is shown in Figure 8 to be more apparent in the PWL forest plot than the FT forest plot. This may be attributed to the taller trees and higher canopy coverage of the PWL forest plot compared to the FT forest plot, as for the same trajectory angle a taller tree will produce a larger downwind footprint.
Figure 9 presents two hemisphere plots which illustrate the correlation between \(C_p\) and interception efficiency at a 0.25 m grid cell resolution over differing azimuth and zenith angles for both the FT and PWL forest plots. These plots demonstrate a strong linear correlation between \(C_p\) and interception efficiency towards the southern portion of the hemisphere, aligning with the average event wind direction. For the PWL forest plot, the upper 97.5th percentile of the \(\rho_p\) values shown in Figure 9, were found between azimuth angles of 167° – 217°. Similarly, for the FT forest plot, the upper 97.5th percentile of \(\rho_p\) was found between azimuth angles of 171°–223°. The zenith angle found to have the highest correlation over this azimuth range was 22° (\(\rho_p\) = 0.7) and 21° (\(\rho_p\) = 0.83) for PWL and FT respectively. The high correlation coefficients found for non-vertical zenith angles for both PWL and FT are hypothesized to result from non-vertical hydrometeor trajectories.
The correlation between \(C_p\) and interception efficiency, resampled to a 5 m grid resolution, was higher when \(C_p\) was adjusted for the observed shift in hydrometeor trajectory (Vector Based), compared to the leaf contact angle measured at a zenith angle of 0° (Nadir) Figure 10. The the zenith angle observed to have the highest \(\rho_p\) in Figure 9 were used to adjust the vector based, \(C_p\) in Figure 10. The stronger association for the vector-based calculation suggests that adjusted \(C_p\) is a useful predictor of interception efficiency before ablation. An ordinary least squares linear regression forced through the origin was fit to the observed data points using the following equation:
\[ \frac{I}{P} = C_p(C_c, \theta_h) \cdot \alpha \tag{9}\]
where \(\alpha\) is an efficiency constant which determines the fraction of snowflakes that contact the \(C_p\) elements and are stored in the canopy (i.e., intercepted) before canopy snow unloading or ablation processes begin.
The Nadir linear regression provided a good overall fit to the observed data with a \(\alpha\) value of 0.95 and 0.99 for the PWL and FT plot respectively (Figure 10). For the PWL plot, the observed points follow the linear relationship up to a \(C_p\) value of around 0.50, above which the increase in interception efficiency levels out. After the Kozak & Kozak (1995) adjustment, a negative R2 value was determined for the PWL plot. Some of the scatter observed in the Nadir model shown in Figure 10 may be explained by grid cells which observed a greater interception efficiency compared to the corresponding \(C_c\) value and can be attributed to the inability of \(C_c\) to represent the increase in interception observed within canopy gaps in Figure 8. Conversely, grid cells where interception efficiency is less than \(C_c\), may be affected by non-vertical trajectory hydrometeors making their way underneath the canopy as observed by the reduced interception efficiency on the windward edges of individual trees in Figure 8. The latter explanation suggests the non-linear relationship observed for the PWL Nadir calculation in Figure 8.
For the vector-based model, the relationship between interception efficiency and \(C_p\) results in R2 values of 0.47 and 0.8 for PWL and FT respectively. The increase in interception efficiency with \(C_p\) follows a reduced slope compared to the Nadir models with \(\alpha\) values of 0.71 and 0.68 for the PWL and FT plots respectively. The reduced slope for the vector-based models may be due to snowflakes that weaved through and/or bounced off branch elements in addition to UAV-lidar measurement uncertainty which may have been slightly affected by unloading and redistribution. These processes would have reduced the fraction of snowfall that was stored in the canopy. Model error statistics are presented in Table 5 for the Nadir and vector-based models and show the vector-based model provided a better prediction of interception efficiency. However, the detailed point clouds required to derive the \(C_p\) values used in this analysis are rarely available and thus more accessible methods to estimate \(C_p\) must be obtained to use Equation 9.
| Plot Name | Canopy Calculation | Model Slope | Mean Bias (-) | MAE (-) | RMS Error (-) | R^2 |
|---|---|---|---|---|---|---|
| FT | Nadir | 0.993 | 0.022 | 0.071 | 0.099 | 0.51 |
| FT | Vector Based | 0.676 | 0.001 | 0.047 | 0.062 | 0.80 |
| PWL | Nadir | 0.949 | 0.048 | 0.113 | 0.146 | NA |
| PWL | Vector Based | 0.707 | 0.019 | 0.078 | 0.095 | 0.47 |
The combined influence of trajectory angle and forest structure on interception
Figure 11 shows that \(C_p\), measured from VoxRS prior to snowfall on March 13th, increases substantially with simulated hydrometeor trajectory angle and corresponding simulated wind speed. The standard deviation in VoxRS measured \(C_p\), illustrated by the shaded area in Figure 11, exhibits the broad range in values for individual grid cells across each forest plot. Despite this large scatter, a systematic increase in the mean \(C_p\) across both forest plots results from a rise in the number of canopy elements for more horizontal portions of the hemisphere, when averaged across each forest plot, over all azimuth angles as shown by the solid lines in the top row of Figure 11. This results in a large rise in VoxRS measured \(C_p\) over relatively common estimated wind speeds. For example, with a wind speed of 1 m s-1 and estimated trajectory angle of 48°, \(C_p\) would increase by 0.31 and 0.28 for the PWL and FT forest plots respectively (Figure 11). This is a fractional increase in the plot \(C_p\) from nadir of 0.61 and 0.95 for PWL and FT respectively. The increase in \(C_p\) from \(C_c\), with increasing trajectory angle is shown on the bottom row of Figure 11 and exhibits a similar relationship for both forest plots FT and PWL until trajectory angles reach approximately 60°. Beyond 60°, the PWL rate of increase slows as the \(C_p\) approaches 1.0, while the FT plot, which has lower \(C_c\), continues to rise until around 75° as a \(C_p\) of 1.0 is approached. \(C_p\) was also quantified across trajectory angles for both PWL and FT on March 14th, post snowfall, and showed a negligible effect from canopy snow load.
A function is proposed here to calculate plot scale leaf contact area, \(C_p\) (-):
\[ C_p = C_c + C_{inc}(\theta_h) \tag{10}\]
where, \(C_{inc}\) is the increase in leaf contact area from \(C_c\) which is a function of \(\theta_h\). To estimate \(C_{inc}\) a non-linear least squares regression using a logistic function forced through the origin was fit to the VoxRS measurements at FT and PWL for simulated hydrometeor trajectory angles (see dashed lines in bottom row of Figure 11). A logistic function was selected to model this relationship, as its shape reflects the slow increase in observed \(C_p\) at near vertical trajectory angles, followed by a rapid increase as snowflakes encounter more canopy area in the middle and lower section of individual trees, and the gradual leveling off of \(C_p\) as it approaches full coverage (value of 1.0). The logistic function used to predict \(C_{inc}\) as a function of \(\theta_h\) is:
\[ C_{inc} = \left( \frac{C^{max}_{inc}}{1 + e ^ {\left(\frac{\theta_0 - \theta_h}{\text{k}}\right)}} - \frac{C^{max}_{inc}}{1 + e ^ {\left(\frac{\theta_0}{\text{k}}\right)}} \right) \tag{11}\]
where \(C^{max}_{inc}\) is the maximum value of \(C_{inc}\), \(\theta_0\) is the x-value of the sigmoid midpoint and \(k\) is the logistic growth rate or steepness of the curve. The coefficients resulting from the non-linear least squares regression fit of Equation 11 to the VoxRS dataset are presented in Table Table 6. Simulated \(C_p\) using Equation 10 is shown in the dashed lines in the top row of Figure 11 and follows closely to the VoxRS-measured mean \(C_p\). Model error statistics shown in Table 7 demonstrate that Equation 11 performed well, with a mean bias and RMSE of 0.001 (-) and 0.0054 (-) respectively for PWL, and -0.0004 (-) and 0.0079 (-) for FT. In contrast, Table 7 reveals that the Hedstrom & Pomeroy (1998) method produced significantly less accurate estimates of \(C_p\), with a mean bias and RMSE of -0.201 (-) and 0.233 (-) respectively for PWL, and -0.260 (-) and 0.324 (-) for FT.
| Plot Names | \(LCA_{max}\) | \(\theta_0\) | \(k\) |
|---|---|---|---|
| PWL | 0.66 | 34.58 | 22.14 |
| FT | 1.18 | 69.13 | 26.98 |
| Model | Plot | Mean Bias (-) | MAE (-) | RMS Error (-) | R^2 |
|---|---|---|---|---|---|
| HP98 | FT | -0.2598 | 0.2598 | 0.3240 | 0.7196 |
| HP98 | PWL | -0.2008 | 0.2010 | 0.2326 | 0.4446 |
| nls | FT | -0.0004 | 0.0067 | 0.0079 | 0.9987 |
| nls | PWL | 0.0010 | 0.0040 | 0.0054 | 0.9990 |
Throughfall Model Performance
The performance of Equations 9, 10, and 11 in estimating event throughfall was assessed against UAV-lidar measurements of throughfall for the March 13–14th snowfall event at the plot scale for both FT and PWL. Required values for the model included the mean hydrometeor terminal velocity and total event snowfall was measured at PWL station, and wind speed was taken as one-third the mean canopy height using the wind speed profile in Figure 7. Additional model inputs include the mean \(C_c\) for each plot which was measured from the VoxRS dataset. An \(\alpha\) value of 0.836 (-) was found through calibration which provided the best fit between observed and simulated interception efficiency at the plot scale for both FT and PWL.
Figure 12 shows the vector-based model, computed using Equation 9 with \(C_p\) adjusted for estimated hydrometeor trajectory angle, closely matches UAV-lidar measurements of throughfall. Observed and modelled values of interception efficiency and \(\Delta SWE_{tf}\) are presented in Table 8 along with corresponding error statistics. Modelled throughfall from the vector-based model was 17 kg m-2 compared to the measured throughfall of 16.6 kg m-2 for PWL. For FT, the modelled throughfall was 21.8 kg m-2, while the measured values where 22.1 kg m-2. The vector-based model shows a lower mean bias of -0.3 kg m-2 for PWL and a negative bias of 0.3 kg m-2 for FT, compared to the larger mean bias of -1.6 kg m-2 for PWL and -0.8 kg m-2 for FT with the nadir-model (calculated using \(C_c\) in place of \(C_p\)). Additionally, the vector-based model significantly reduced error in predicted throughfall, from -9.4% with the nadir-model to -1.8% for PWL. A smaller improvement was observed for FT, with the error in predicted throughfall declining from -3.6% with the nadir-model to -1.4% with the vector-based model.
| Plot | Model Type | Value Name | Units | Obs. Value | Mod. Value | Mean Bias | Perc. Error |
|---|---|---|---|---|---|---|---|
| FT | VB-model | I/P | - | 0.23 | 0.24 | -0.01 | -4.67 |
| FT | Nadir-model | I/P | - | 0.23 | 0.20 | 0.03 | 12.10 |
| FT | VB-model | TF | kg m⁻² | 22.12 | 21.82 | 0.31 | 1.38 |
| FT | Nadir-model | TF | kg m⁻² | 22.12 | 22.91 | -0.79 | -3.58 |
| PWL | VB-model | I/P | - | 0.42 | 0.41 | 0.01 | 2.54 |
| PWL | Nadir-model | I/P | - | 0.42 | 0.37 | 0.05 | 12.95 |
| PWL | VB-model | TF | kg m⁻² | 16.64 | 16.95 | -0.31 | -1.84 |
| PWL | Nadir-model | TF | kg m⁻² | 16.64 | 18.20 | -1.56 | -9.35 |
Discussion
The point scale observations presented in Figure 6 show air temperature had little influence on interception efficiency. This differs from existing studies which suggested either a positive (Storck et al., 2002) or negative (Hedstrom & Pomeroy, 1998) relationship. An increase in initial interception efficiency before unloading was observed with increasing wind speed at two locations which were sheltered from the predominant wind direction (Figure 6, B). This is attributed to an associated increase in \(C_p\) with wind speed. These results are consistent with observations by Schmidt & Troendle (1989) who observed a slight increase in snowfall interception with increasing wind speeds up to 6 m s-1 and studies of rainfall interception by Herwitz & Slye (1995) and Van Stan et al. (2011).
Compared to the influence of wind speed, interception efficiency showed a smaller sensitivity to canopy snow load at the point scale (Figure 5). The slight increase in interception efficiency for smaller canopy snow loads and decline for larger canopy snow loads is attributed to the influence of canopy snow load on \(C_p\) (Figure 6, C). While small, this effect is consistent with the theory proposed by Satterlund & Haupt (1967) that interception efficiency increases as the canopy fills with snow bridging gaps in the canopy increasing, while later declining due to branch bending and decreased canopy coverage. However, the observations presented in Figure 6 and Figure 3, differ from the Satterlund & Haupt (1967), Hedstrom & Pomeroy (1998), Storck et al. (2002) and Moeser et al. (2015) theories, as canopy snow load increased linearly with snowfalls up to 45 kg m-2 with no evidence of approaching a maximum canopy snow load. The strong exponential decline in interception efficiency observed with increasing event snowfall in Satterlund & Haupt (1967), Hedstrom & Pomeroy (1998), Storck et al. (2002) and Moeser et al. (2015) may have been a result of increased unloading rates as branches bent under heavy snow loads and hence mixed ablation and interception processes to varying degrees. The low sensitivity of interception efficiency with canopy snow load found in this study may be attributed to several factors: a reduced inclusion of ablation processes in the interception efficiency measurements, limited influence of canopy snow load on \(C_p\) at this study site, and/or the compensatory effects outlined by Satterlund & Haupt (1967).
Staines & Pomeroy (2023) showed a slight increase in VoxRS \(C_p\) between snow-off and snow-on conditions. However, the increase in \(C_p\) resulting from snow load in Staines & Pomeroy (2023) was small compared to the substantial rise in \(C_p\) due to trajectory angle presented in their study and as shown in Figure 11. Both findings from Staines & Pomeroy (2023) corroborate the results reported in this study. Further evidence in support of the relatively small influence of canopy snow load on \(C_p\), is provided by Lundquist et al. (2021) who reported improved simulation of subcanopy snow accumulation without the use of a maximum canopy snow load, when linked with a comprehensive canopy snow ablation routine. Lehtonen et al. (2016) noted that in northern Finland heavy canopy snow loads have been observed to continue increasing until stem breakage, under conditions favourable for the formation of significant rime-ice accretion and limited ablation, thus reducing \(C_p\). Models are available to predict the accretion of ice on tree canopies (e.g., Nock et al., 2016) however, further research is required to understand the canopy snow load required to cause stem breakage across different tree species and canopy loads.
These findings on the limited influence of air temperature and canopy snow load on initial interception challenge the theoretical basis of many existing snow interception parameterizations (Hedstrom & Pomeroy, 1998; Moeser et al., 2015; Satterlund & Haupt, 1967; Storck et al., 2002). To address this a new snow interception parameterization, Equation 9, is presented which calculates interception efficiency as a function of \(C_p\) and \(\alpha\). This new parameterization allows for canopy snow loading processes to be isolated from canopy snow ablation processes and is consistent with current rainfall interception theory (Valante et al., 1997). Equation 9 differs only slightly from the original Hedstrom & Pomeroy (1998) parameterization (see Equation 6 in Hedstrom & Pomeroy (1998)), in that it does not calculate interception efficiency as a function of canopy snow load and from the Storck et al. (2002) parameterization who proposed interception efficiency to be constant over time and space. The theoretical basis of the \(\alpha\) parameter in Equation 9 is that the association between \(C_p\) and interception efficiency, as shown in Figure 10, does not follow a 1:1 line as falling snow hydrometeors may bounce off the canopy elements. Further research is needed to explore how processes such as the increased cohesion and adhesion of snowfall to the canopy at warm temperatures, as observed by Kobayashi (1987), Pfister & Schneebeli (1999), Storck et al. (2002), as well as hydrometeor velocity, particle size, and shape suggested by (Katsushima et al., 2023), may influence the \(\alpha\) parameter, although these effects were not observed in this study.
Measurements of interception efficiency and canopy structure, as shown in Figure 8, align with the theory proposed by Hedstrom & Pomeroy (1998) which suggests reduced throughfall on the lee side of individual trees. However, an existing method proposed in Hedstrom & Pomeroy (1998) to scale canopy coverage with wind speed failed to reproduce the observations presented in Figure 11. A new method is proposed which uses a logistic function to calculate plot-wide \(C_{inc}\) as a function of \(\theta_h\) and \(C_c\). Significant scatter in VoxRS measured \(C_p\) across the two forest plots, illustrated by the high standard deviation in Figure 11, resulted from directional (azimuth) and spatial differences in canopy structure. This large scatter suggests the observed relationships in Figure 11 are only applicable at the forest stand scale where the sub metre variability in \(C_p\) averages out. At the point scale, the mixed canopy SCL which is open to the prevailing wind direction (Figure 2), and did not follow this relationship and led to an increase in throughfall with increasing wind speed (Figure 5 & Figure 6). However, Figure 11 shows that at the plot scale, while some grid cells with an opening towards a certain azimuth decreasing \(C_p\) with increasing \(\theta_h\), there is a greater number of grid cells which have closed canopy, driving the plot wide increase in \(C_p\) with \(\theta_h\). Thus at the plot scale, Equation 11, which uses trajectory angle alone, could be used to determine \(C_{inc}\) and thus \(C_p\) even for the discontinuous canopies in the FT and PWL forest plots. However, Equation 11 would not be applicable to areas that have large continuous gap fractions (e.g., large forested clear cuts). Further work is required to refine the relationship proposed in Equation 11 across a range of tree species and densities. It was found that the mean hydrometeor trajectory angle over a snowfall event, required for Equation 11, could be predicted by using the observed hydrometeor fall velocity and a mean horizontal wind speed selected at one-third of the canopy height above the ground. A wind speed at one-third the mean canopy height is hypothesized to be important for canopy snow accumulation as a large fraction of the horizontal cross-sectional area is at this height for most needleleaf canopies. Katsushima et al. (2023), also proposed the wind speed at one-third the canopy height for modelling unloading of canopy snow as it corresponds to the centre of gravity when the horizontal projection of the canopy is assumed to be a triangle. However, there is uncertainty in the transferability of one-third canopy height observed here to other environments due to differing tree structures and tree species such as those with a larger trunk space or have more of their canopy contact area at higher heights above the ground (i.e., some deciduous canopies). Moreover, Equation 4 assumes a linear hydrometeor trajectory, and does not consider non-linear patterns such as wind flow directions around tree elements, turbulent flow, or differences in wind speed with height.
Although the improvement in performance of the vector-based model over the Nadir model was relatively small, the vector-based model is preferred due to its overall lower error compared to the UAV-lidar measurements and better representation of physical processes. While the vector-based model acts to increase interception efficiency with wind speed, several studies have shown that canopy snow ablation increases as a result of wind induced unloading (Bartlett & Verseghy, 2015; Betts & Ball, 1997; Lumbrazo et al., 2022; Roesch et al., 2001; Wheeler, 1987). While these studies have been used to develop parameterizations for wind induced unloading, they were not based on direct measurements of canopy snow unloading. Further research and direct in-situ observations of canopy snow ablation are required to refine ablation parameterization.
Conclusions
New observations of initial snow interception, collected over a wide range of meteorological conditions and canopy structures suggest forest structure is the primary factor governing subcanopy snow accumulation. At the point scale, high-temporal resolution measurements revealed no evidence of a maximum canopy snow load, even for event snowfalls up to 45 kg m-2, nor was there any indication of air temperature influencing the cohesion and adhesion of snowfall to the canopy or branch bending reducing canopy coverage. Instead, wind speed was found to influence interception efficiency by changing the hydrometeor trajectory angle, which can lead to a substantial increase in apparent forest structure or snow-leaf contact area.
At the forest plot scale, UAV-lidar measurements of throughfall collected over a wind-driven snowfall event confirmed the results observed at the point-scale and showed leaf contact area was the main factor governing the interception efficiency at a particular site. The leaf contact area, which accounts for the change in canopy structure with trajectory angle, proved to be a better predictor of interception efficiency compared to nadir-calculated canopy coverage. When averaged across each forest plot, leaf contact area was shown to be highly sensitive to hydrometeor trajectory angle, increasing by 61–95% for trajectory angles associated with a 1 m s-1 wind speed. An existing theoretical relationship failed to adequately represent the VoxRS-measured increase in leaf contact area with simulated trajectory angles. As a result, a new relationship is proposed, which demonstrated good performance at this study site.
The weak association between air temperature and canopy snow load with interception efficiency, as presented here and in other recent studies, coupled with the considerable influence of wind speed on leaf contact area, highlights the need for a new snow interception parameterization. A new parameterization is proposed that calculates initial interception as a function of snowfall and leaf contact area. This parameterization is consistent with many rainfall interception studies, which also separate canopy loading and ablation processes and calculate interception as a function of canopy coverage. Additionally, a second equation is proposed to estimate leaf contact area as a function of hydrometeor trajectory angle and nadir canopy coverage. This updated snow interception parameterization showed good performance in the subalpine forest in this study, but further validation should be conducted in a range of climates, forest species, and canopy structures.
Acknowledgments
We wish to acknowledge financial support from the University of Saskatchewan Dean’s Scholarship, Natural Sciences and Engineering Research Council of Canada, Global Water Futures Programme, Alberta Innovates, Canada Foundation for Innovation, and the Canada Research Chairs Programme. We thank Madison Harasyn, Hannah Koslowsky, Kieran Lehan, Lindsey Langs and Fortress Mountain Resort for their help in the field. Madison Harasyn, Alistair Wallace, and Rob White contributed to developing the UAV-lidar processing workflow.
Data Availability
The data that support the findings in this study will be made publicly available upon publication.
References
Axelsson, P. (2000). Peter Axelsson 110. International Archives of Photogrammetry and Remote Sensing, 33(4), 110–117. https://www.isprs.org/proceedings/XXXIII/congress/part4/111_XXXIII-part4.pdf
Bartlett, P. A., & Verseghy, D. L. (2015). Modified treatment of intercepted snow improves the simulated forest albedo in the Canadian Land Surface Scheme. Hydrological Processes, 29(14), 3208–3226. https://doi.org/10.1002/HYP.10431
BayesMap Solutions. (2024). BayesStripAlign. https://bayesmap.com/products/bayesstripalign/
Betts, A. K., & Ball, J. H. (1997). Albedo over the boreal forest. Journal of Geophysical Research: Atmospheres, 102(D24), 28901–28909. https://doi.org/https://doi.org/10.1029/96JD03876
Chianucci, F., & Macek, M. (2023). hemispheR: an R package for fisheye canopy image analysis. Agricultural and Forest Meteorology. https://doi.org/10.1016/j.agrformet.2023.109470
Cionco, R. M. (1965). A Mathematical Model for Air Flow in a Vegetative Canopy. Journal of Applied Meteorology (1962), 4(4), 517–522. https://doi.org/https://doi.org/10.1175/1520-0450(1965)004<0517:AMMFAF>2.0.CO;2
Clark, M. P., Wood, A., Nijssen, B., Bennett, A., Knoben, W., & Lumbrazo, C. (2020). SUMMA v3.0.3. Zenodo. https://doi.org/10.5281/zenodo.4558054
Deems, J. S., Painter, T. H., & Finnegan, D. C. (2013). Lidar measurement of snow depth: A review. Journal of Glaciology, 59(215), 467–479. https://doi.org/10.3189/2013JoG12J154
Ellis, C. R., Pomeroy, J. W., & Link, T. E. (2013). Modeling increases in snowmelt yield and desynchronization resulting from forest gap-thinning treatments in a northern mountain headwater basin. Water Resources Research, 49(2), 936–949. https://doi.org/10.1002/wrcr.20089
Essery, R., Pomeroy, J. W., Parviainen, J., & Storck, P. (2003). Sublimation of snow from coniferous forests in a climate model. Journal of Climate, 16(11), 1855–1864. https://doi.org/10.1175/1520-0442(2003)016<1855:SOSFCF>2.0.CO;2
Floyd, W. C. (2012). Snowmelt energy flux recovery during rain-on-snow in regenerating forests (p. 180) [PhD Thesis, University of British Columbia]. https://doi.org/https://dx.doi.org/10.14288/1.0073024
Fryer, B. Y. G. I., Johnson, E. A., Fryer, G. I., & Johnson, E. A. (1988). Reconstructing Fire Behaviour and Effects in a Subalpine Forest. The Journal of Applied Ecology, 25(3), 1063–1072. https://doi.org/https://doi.org/10.2307/2403766
Gelfan, A. N., Pomeroy, J. W., & Kuchment, L. S. (2004). Modeling forest cover influences on snow accumulation, sublimation, and melt. Journal of Hydrometeorology, 5(5), 785–803. https://doi.org/10.1175/1525-7541(2004)005<0785:MFCIOS>2.0.CO;2
Golding, D. L., & Swanson, R. H. (1978). Snow accumulation and melt in small forest openings in Alberta. Canadian Journal of Forest Research, 8(4), 380–388. https://doi.org/https://doi.org/10.1139/x78-057
Harder, P., & Pomeroy, J. W. (2013). Estimating precipitation phase using a psychrometric energy balance method. Hydrological Processes, 27(13), 1901–1914. https://doi.org/https://doi.org/10.1002/hyp.9799
Harder, P., Pomeroy, J. W., & Helgason, W. D. (2020). Improving sub-canopy snow depth mapping with unmanned aerial vehicles: Lidar versus structure-from-motion techniques. Cryosphere, 14(6), 1919–1935. https://doi.org/10.5194/tc-14-1919-2020
Harpold, A. A., Krogh, S. A., Kohler, M., Eckberg, D., Greenberg, J., Sterle, G., & Broxton, P. D. (2020). Increasing the efficacy of forest thinning for snow using high-resolution modeling: A proof of concept in the Lake Tahoe Basin, California, USA. Ecohydrology, 13(4). https://doi.org/10.1002/eco.2203
Hedstrom, N. R., & Pomeroy, J. W. (1998). Measurements and modelling of snow interception in the boreal forest. Hydrological Processes, 12(10-11), 1611–1625. https://doi.org/10.1002/(SICI)1099-1085(199808/09)12:10/11<1611::AID-HYP684>3.0.CO;2-4
Herwitz, S. R., & Slye, R. E. (1995). Three-dimensional modeling of canopy tree interception of wind-driven rainfall. Journal of Hydrology, 168(1-4), 205–226. https://doi.org/10.1016/0022-1694(94)02643-P
Hijmans, R. J. (2024). terra: Spatial Data Analysis. https://cran.r-project.org/package=terra
Katsushima, T., Kato, A., Aiura, H., Nanko, K., Suzuki, S., Takeuchi, Y., & Murakami, S. (2023). Modelling of snow interception on a Japanese cedar canopy based on weighing tree experiment in a warm winter region. Hydrological Processes, 37(6), 1–16. https://doi.org/10.1002/hyp.14922
Kim, E., Gatebe, C., Hall, D., Newlin, J., Misakonis, A., Elder, K., Marshall, H. P., Hiemstra, C., Brucker, L., De Marco, E., Crawford, C., Kang, D. H., & Entin, J. (2017). NASA’s snowex campaign: Observing seasonal snow in a forested environment. 2017 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), 1388–1390. https://doi.org/10.1109/IGARSS.2017.8127222
Kobayashi, D. (1987). Snow accumulation on a narrow board. Cold Regions Science and Technology, 13(3), 239–245. https://doi.org/https://doi.org/10.1016/0165-232X(87)90005-X
Kozak, A., & Kozak, R. A. (1995). Notes on regression through the origin. Forestry Chronicle, 71(3), 326–330. https://doi.org/https://doi.org/10.5558/tfc71326-3
Langs, L. E., Petrone, R. M., & Pomeroy, J. W. (2020). A \(\delta\)18O and \(\delta\)2H stable water isotope analysis of subalpine forest water sources under seasonal and hydrological stress in the Canadian Rocky Mountains. Hydrological Processes, 34(26), 5642–5658. https://doi.org/10.1002/hyp.13986
LAStools. (2024). Efficient LiDAR Processing Software (version 240220, academic). http://rapidlasso.com/LAStools
Lehtonen, I., Kämäraïnen, M., Gregow, H., Venälaïnen, A., & Peltola, H. (2016). Heavy snow loads in Finnish forests respond regionally asymmetrically to projected climate change. Natural Hazards and Earth System Sciences, 16(10), 2259–2271. https://doi.org/10.5194/nhess-16-2259-2016
Lumbrazo, C., Bennett, A., Currier, W. R., Nijssen, B., & Lundquist, J. (2022). Evaluating Multiple Canopy-Snow Unloading Parameterizations in SUMMA With Time-Lapse Photography Characterized by Citizen Scientists. Water Resources Research, 58(6), 1–22. https://doi.org/10.1029/2021WR030852
Lundberg, A., & Hallidin, S. (1994). Evaporation of intercepted snow: Analysis of governing factors. Water Resources Research, 30(9), 2587–2598. https://doi.org/10.1029/94WR00873
Lundquist, J. D., Dickerson-Lange, S., Gutmann, E., Jonas, T., Lumbrazo, C., & Reynolds, D. (2021). Snow interception modelling: Isolated observations have led to many land surface models lacking appropriate temperature sensitivities. Hydrological Processes, 35(7), 1–20. https://doi.org/10.1002/hyp.14274
MacDonald, J. P. J. (2010). Unloading of intercepted snow in conifer forests (Master of Science August, University of Saskatchewan; pp. 0–93). http://hdl.handle.net/10388/etd-09302010-145706
Mahat, V., & Tarboton, D. G. (2014). Representation of canopy snow interception, unloading and melt in a parsimonious snowmelt model. Hydrological Processes, 28(26), 6320–6336. https://doi.org/10.1002/hyp.10116
Moeser, D., Stähli, M., & Jonas, T. (2015). Improved snow interception modeling using canopy parameters derived from airborne LiDAR data. Water Resources Research, 51(7), 5041–5059. https://doi.org/https://doi.org/10.1002/2014WR016724
Musselman, K. N., Pomeroy, J. W., & Link, T. E. (2015). Variability in shortwave irradiance caused by forest gaps: Measurements, modelling, and implications for snow energetics. Agricultural and Forest Meteorology, 207, 69–82. https://doi.org/https://doi.org/10.1016/j.agrformet.2015.03.014
Natural Resources Canada. (2024). Precise point positioning. https://webapp.csrs-scrs.nrcan-rncan.gc.ca/geod/tools-outils/ppp.php
Nock, C. A., Lecigne, B., Taugourdeau, O., Greene, D. F., Dauzat, J., Delagrange, S., & Messier, C. (2016). Linking ice accretion and crown structure: towards a model of the effect of freezing rain on tree canopies. Annals of Botany, 117(7), 1163–1173. https://doi.org/https://doi.org/10.1093/aob/mcw059
Pfister, R., & Schneebeli, M. (1999). Snow accumulation on boards of different sizes and shapes. Hydrological Processes, 13(14-15), 2345–2355. https://doi.org/https://doi.org/10.1002/(SICI)1099-1085(199910)13:14/15<2345::AID-HYP873>3.0.CO;2-N
Pomeroy, J. W., Brown, T., Fang, X., Shook, K. R., Pradhananga, D., Armstrong, R., Harder, P., Marsh, C., Costa, D., Krogh, S. A., Aubry-wake, C., Annand, H., Lawford, P., He, Z., Kompanizare, M., & Moreno, J. I. L. (2022). The cold regions hydrological modelling platform for hydrological diagnosis and prediction based on process understanding. Journal of Hydrology, 615(128711), 1–25. https://doi.org/10.1016/j.jhydrol.2022.128711
Pomeroy, J. W., Fang, X., & Ellis, C. R. (2012). Sensitivity of snowmelt hydrology in Marmot Creek, Alberta, to forest cover disturbance. Hydrological Processes, 26(12), 1891–1904. https://doi.org/10.1002/hyp.9248
Pomeroy, J. W., Marsh, P., & Gray, D. M. (1997). Application of a distributed blowing snow model to the arctic. Hydrological Processes, 11(11), 1451–1464. https://doi.org/10.1002/(sici)1099-1085(199709)11:11<1451::aid-hyp449>3.0.co;2-q
Pomeroy, J. W., Parviainen, J., Hedstrom, N., & Gray, D. M. (1998). Coupled modelling of forest snow interception and sublimation. Hydrological Processes, 12(15), 2317–2337. https://doi.org/10.1002/(SICI)1099-1085(199812)12:15<2317::AID-HYP799>3.0.CO;2-X
Pomeroy, J. W., & Schmidt, R. A. (1993). The Use of Fractal Geometry in Modelling Intercepted Snow Accumulation and Sublimation. Eastern Snow Conference, 50, 231–239.
R Core Team. (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.r-project.org/
Rittger, K., Raleigh, M. S., Dozier, J., Hill, A. F., Lutz, J. A., & Painter, T. H. (2020). Canopy Adjustment and Improved Cloud Detection for Remotely Sensed Snow Cover Mapping. Water Resources Research, 56(6), n/a. https://doi.org/10.1029/2019WR024914
Roesch, A., Wild, M., Gilgen, H., & Ohmura, A. (2001). A new snow cover fraction parameterization for the ECHAM4 GCM. Climate Dynamics, 17(12), 933–946. https://doi.org/10.1007/s003820100153
Safa, H., Krogh, S. A., Greenberg, J., Kostadinov, T. S., & Harpold, A. A. (2021). Unraveling the Controls on Snow Disappearance in Montane Conifer Forests Using Multi-Site Lidar. Water Resources Research, 57(12), 1–20. https://doi.org/10.1029/2020WR027522
Satterlund, D. R., & Haupt, H. F. (1967). Snow catch by Conifer Crowns. Water Resources Research, 3(4), 1035–1039. https://doi.org/https://doi.org/10.1029/WR003i004p01035
Schmidt, R. A., & Troendle, C. A. (1989). Snowfall into a forest and clearing. Journal of Hydrology, 110(3-4), 335–348. https://doi.org/10.1016/0022-1694(89)90196-0
Smith, C. D. (2007). Correcting the wind bias in snowfall measurements made with a Geonor T-200B precipitation gauge and alter wind shield. 87th AMS Annual Meeting.
Staines, J., & Pomeroy, J. W. (2023). Influence of forest canopy structure and wind flow on patterns of sub-canopy snow accumulation in montane needleleaf forests. Hydrological Processes, 37(10), 1–19. https://doi.org/10.1002/hyp.15005
Storck, P., Lettenmaier, D. P., & Bolton, S. M. (2002). Measurement of snow interception and canopy effects on snow accumulation and melt in a mountainous maritime climate, Oregon, United States. Water Resources Research, 38(11), 1–16. https://doi.org/10.1029/2002wr001281
Troendle, C. A. (1983). The Potential for Water Yield Augmentation From Forest Management in the Rocky Mountain Region. Journal of the American Water Resources Association, 19(3), 359–373. https://doi.org/10.1111/j.1752-1688.1983.tb04593.x
Valante, F., David, J. S., & Gash, J. H. C. (1997). Modelling interception loss for two sparse eucalypt and pine forests in central Portugal using reformulated Rutter and Gash analytical models. Journal of Hydrology, 190(1-2), 141–162. https://doi.org/10.1016/S0022-1694(96)03066-1
Van Stan, J. T., Siegert, C. M., Levia, D. F., & Scheick, C. E. (2011). Effects of wind-driven rainfall on stemflow generation between codominant tree species with differing crown characteristics. Agricultural and Forest Meteorology, 151(9), 1277–1286. https://doi.org/10.1016/j.agrformet.2011.05.008
Varhola, A., Coops, N. C., Weiler, M., & Moore, R. D. (2010). Forest canopy effects on snow accumulation and ablation: An integrative review of empirical results. Journal of Hydrology, 392(3-4), 219–233. https://doi.org/10.1016/j.jhydrol.2010.08.009
Vionnet, V., Mortimer, C., Brady, M., Arnal, L., & Brown, R. (2021). Canadian historical Snow Water Equivalent dataset (CanSWE, 1928–2020). Earth System Science Data, 13(9), 4603–4619. https://doi.org/https://doi.org/10.5194/essd-13-4603-2021
Wheeler, K. (1987). Interception and redistribution of snow in a subalpine forest on a storm-by-storm basis. Western Snow Conference. http://sites/westernsnowconference.org/PDFs/1987Wheeler.pdf
Supporting Information
Detailed Description of UAV-Lidar Methodology
The REIGL miniVUX-2 laser operates at a near infrared wavelength with a laser beam footprint of 0.160 m x 0.05 mm (at 100 m above ground). The accuracy and precision of the miniVUX-2 is described by REIGL for a lab environment of 0.015 m and 0.01 m respectively (at 50 m above ground). The miniVUX-2 was configured with a laser pulse repetition rate of 200 kHz, field of view of 360°, scan speed of 31.09 revolutions s-1 and an angular step width of 0.0558°, resulting in an expected an average point cloud density of 107 returns m-2 for each flight path.
Georeferenced point clouds with x, y, and z coordinates for each laser return were generated following methods outlined by Harder et al. (2020) and Staines & Pomeroy (2023) to reconcile survey lidar, IMU and GNSS data. A ground-based GNSS system was positioned on a permanent monument during each survey and underwent precise point positioning (PPP) correction by Natural Resources Canada (2024). Differential GNSS correction of the UAV trajectory was conducted using the ground-based PPP GNSS observations and the POSPac UAV software. The UAV-lidar point clouds were then transformed from a sensor referenced coordinate system to a georeferenced coordinate system (EPSG:32611 - WGS 84 / UTM zone 11N) using the RIEGL Riprocess Software. A vertical offset of up to 6 cm between UAV-lidar flight lines was observed in the resulting point clouds on March 13th and 14th, 2024 and was attributed to IMU position drift. This offset between flight lines was corrected using the BayesStripAlign software v2.24 (BayesMap Solutions, 2024), which reduces relative and absolute uncertainties in the vertical elevation of the point cloud using the ground control points (GCP) collected across the study site using a differential GNSS rover.
Quality control, ground classification and calculation of the change in between two UAV-lidar point clouds was conducted using the LAStools software package (LAStools, 2024). The ground classification was conducted using the “lasground_new” function (LAStools, 2024) for both the pre and post snowfall event point clouds, with a step size set to 2 m and 8 substeps (ultra_fine setting). The offset and spike options were set to remove points that are more than 0.1 m above or below the initial ground surface estimate surface which “lasground_new” fits to the last returns. This function is based on an algorithm outlined by Axelsson (2000), describing the process of making the initial ground surface element.
The change in elevation between the two UAV-lidar surveys was interpreted as the increase in snow accumulation, \(\Delta HS\) over the snowfall event. This change was calculated using a point-to-grid subtraction method, using the “lasheight” function from the LAStools (2024) software, as in Deems et al. (2013) and Staines & Pomeroy (2023). The pre snowfall event point cloud from “lasground_new” by “lasheight” to construct a “ground” TIN. Subsequently, the height of each post snowfall event point above the ground TIN, resulting in a point cloud representing \(\Delta HS\). This point cloud was then converted into a raster of \(\Delta HS\) with a grid cell resolution of 5 x 5 cm using the “las2dem” function. Further quality control and resampling of the 5 cm raster of \(\Delta HS\) was conducted using the ‘Terra’ R package (Hijmans, 2024). Areas that were disturbed over the snowfall event during the in-situ snow survey and values that exceeded the .999th quantile were removed. To help remove any remaining noise a 25 cm \(\Delta HS\) raster was generated by computing the median of the 5 cm \(\Delta HS\) values within each 25 cm grid cell.
A comparison of UAV-liar and in-situ snow survey measurements over the March 13–14th snowfall event and associated error metrics are shown in Figure 13.
Linear Regression Models Through the Origin
Kozak & Kozak (1995) noted, the default R2 value provided for least squares models forced through the origin by many statistical packages can be misleading. Therefore, these R2 values were adjusted using Equation 10 in Kozak & Kozak (1995) and two statistical tests as described by Kozak & Kozak (1995) were used to verify whether a no-intercept model (forced through the origin) was appropriate for this data compared to a with-intercept model. The first test evaluated if the intercept of the with-intercept was significantly different from zero using p-value provided by the ‘summary’ function from the ‘stats’ package in R (R Core Team, 2024). The second test examined if there was a significant difference between the no-intercept and with-intercept models by testing if the residual sum of squares was different between the no-intercept and full model, assessed via Equation 15 in Kozak & Kozak (1995). If the first test indicated a significant difference, and the second did not, the no-intercept model could be deemed statistically justified (Kozak & Kozak, 1995).